Writing a User Material for

Commercial Finite Element Software

Lecture adapted from Abaqus lecture on writing UMATs

Overview

• Motivation
• Steps Required for Writing a User Material
• Abaqus UMAT Interface
• Linear Elastic Example

Overview, contd.

Many commercial finite element codes provide an interface that allows users to implement general material constitutive equations. A user material model:

• is used to define the mechanical constitutive behavior of a material;

• is called at integration points of elements for which the material definition includes a user-defined material behavior;

• can use solution-dependent state variables;

• updates the stresses and solution-dependent state variables to their values at the end of the increment for which it is called;

• for most implicit solvers must provide the material Jacobian matrix, $\partial\Delta\pmb{\sigma}/\partial\Delta\pmb{\epsilon}$;

Use the user material interface only when none of the existing material models included in the code's material library accurately represents the behavior of the material to be modeled.

Steps Required in Writing a User Material

Define the Constitutive Response

• Outside the scope of this lecture

• Requires one of the following:

  • Explicit definition of stress (Cauchy stress for large-strain applications)
  • Definition of the stress rate only (in corotational framework)

• It is also likely to require:

  • Definition of dependence on time, temperature, or field variables
  • Definition of internal state variables, either explicitly or in rate form

Transform Constitutive Equations to Incremental Form

Use a suitable integration procedure (this is the hard part!):

• Forward Euler (explicit integration)

  Simple to implement, but has a stability limit

  $$\Delta\epsilon < \Delta\epsilon_{\rm stab}$$

  where $\Delta\epsilon_{\rm stab}$ is usually less than the elastic strain magnitude.

• Backward Euler (implicit integration)

• Midpoint method

Considerations

• For explicit integration the time increment must be controlled.
• For implicit or midpoint integration the algorithm is more complicated and often requires local iteration. However, there is usually no stability limit.
• An incremental expression for the internal state variables must also be obtained.

Calculation of the (consistent) Jacobian

• Required only for implicit solvers

• For small-deformation problems (e.g., linear elasticity) or large-deformation problems with small volume changes (e.g., metal plasticity), the consistent Jacobian is

  $$ \mathbb{C} = \frac{\partial\Delta\pmb{\sigma}}{\partial\Delta\pmb{\epsilon}}$$

  where $\Delta\sigma$ is the increment in stress and $\Delta\epsilon$ is the increment in strain.

• $\pmb{\sigma}$ is most always associated with the Cauchy stress and, in finite-strain problems, $\pmb{\epsilon}$ is an approximation to the logarithmic strain.
• This matrix may be nonsymmetric as a result of the constitutive equation or integration procedure.
• The Jacobian is often approximated, resulting in a loss of quadratic convergence.
• It is easily calculated for forward integration methods (usually the elasticity matrix).

Calculation of the (consistent) Jacobian, contd.

• If large deformations with large volume changes are considered (e.g., pressure-dependent plasticity), the exact form of the consistent Jacobian

  $$ \mathbb{C} = \frac{1}{J}\frac{\partial\Delta\left(J\pmb{\sigma}\right)}{\partial\Delta\pmb{\epsilon}}$$

  should be used to ensure rapid convergence.

• For hyperelastic constitutive equations, in which $\pmb{\sigma}$ is a proper function of the deformation, the consistent Jacobian is defined by:

  $$\delta\left(J\pmb{\sigma}\right) = J\mathbb{C}{:}\delta\pmb{d}$$

Commercial Code Interfaces

Many commercial codes provide one or more user material interfaces

• Abaqus/Standard: UMAT (Fortran)
• Abaqus/Explicit: VUMAT (Fortran)
• LSDYNA: UMAT (Fortran)
• Ansys: USERMAT (Fortran)
• FEBio: Subclass of base material (C++)

Abaqus UMAT Interface

• Used in Abaqus/Standard
• Follows Fortran coding conventions
• User must make sure that all variables are defined and initialized properly.
• Use Abaqus utility routines as required.
• Assign enough storage space for state variables with the ∗DEPVAR option.
• Abaqus UMATs can be run directly in Matmodlab

Verification

• Verify the subroutine response using Matmodlab

• In Abaqus, run tests with all displacements prescribed to verify the integration algorithm for stresses and state variables.

  Suggested tests include:

  • Uniaxial
  • Uniaxial in oblique directions
  • Uniaxial with finite rotation
  • Finite shear

• Run similar tests with load prescribed to verify the accuracy of the Jacobian.

  • In Matmodlab, use the sqa_stiff option to compare stiffness to that computed numerically
• Compare test results with analytical solutions.

Abaqus UMAT Interface

• See the Abaqus documentation for input file requirements
• Abaqus documents that UMATs must be written in Fortran 77, but can be implented in free form Fortran
• If free form Fortran is used, the Fortran command in the Abaqus environment file must be modified
• The UMAT Fortran 77 subroutine header is

SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, SCD, RPL,
      1 DDSDDT, DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,
      2 PREDEF,DPRED, CMNAME, NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS, 
      3 COORDS, DROT, PNEWDT, CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER, 
      4 KSPT, KSTEP, KINC)

       INCLUDE 'ABA_PARAM.INC'
       CHARACTER*8 CMNAME
       DIMENSION STRESS(NTENS), STATEV(NSTATV), DDSDDE(NTENS, NTENS),
      1 DDSDDT(NTENS), DRPLDE(NTENS), STRAN(NTENS), DSTRAN(NTENS),
      2 PREDEF(1), DPRED(1), PROPS(NPROPS), COORDS(3), DROT(3, 3),
      3 DFGRD0(3, 3), DFGRD1(3, 3)

• The include statement sets the proper precision for floating point variables
• If implicit none is used, make sure the precision set matches that of Abaqus

UMAT Variables

The following quantities are available in UMAT:

• Stress, strain, and SDVs at the start of the increment
• Strain increment, rotation increment, and deformation gradient at the start and end of the increment
• Total and incremental values of time, temperature, and user-defined field variables
• Material constants, material point position, and a characteristic element length
• Element, integration point, and composite layer number (for shells and layered solids)
• Current step and increment numbers

The following quantities must be defined:

• Stress, SDVs, and material Jacobian

The following variables may be defined:

• Strain energy, plastic dissipation, and “creep” dissipation
• Suggested new (reduced) time increment

A complete description of all parameters is given in the Abaqus user subroutines guide

UMAT Utilities

• SINV will return the first and second invariants of a tensor.
• SPRINC will return the principal values of a tensor.
• SPRIND will return the principal values and directions of a tensor.
• ROTSIG will rotate a tensor with an orientation matrix.
• XIT will terminate an analysis and close all files associated with the analysis properly.
• STD_ABQERR sends messages to Abaqus to write to output files

For more details, see the Abaqus user subroutines guide

UMAT Conventions

Sresses and strains are stored as vectors

• For plane stress elements: $\sigma_{11}, \sigma_{22}, \sigma_{12}$
• For (generalized) plane strain and axisymmetric elements: $\sigma_{11}, \sigma_{22}, \sigma_{33}, \sigma_{12}$
• For three-dimensional elements: $\sigma_{11}, \sigma_{22}, \sigma_{33}, \sigma_{12}, \sigma_{13}, \sigma_{23}$

The shear strain is stored as engineering shear strain

The deformation gradient, $F_{ij}$, is always stored as a three-dimensional matrix

UMAT Formulation Aspects

• For geometrically nonlinear analysis the strain increment and incremental rotation passed into the routine are based on the Hughes-Winget formulae.

  • Linearized strain and rotation increments are calculated in the mid-increment configuration.
  • Approximations are made, particularly if rotation increments are large: more accurate measures can be obtained from the deformation gradient if desired.

• The user must define the Cauchy stress: when this stress reappears during the next increment, it will have been rotated with the incremental rotation, DROT, passed into the subroutine.

  • The stress tensor can be rotated back using the utility routine ROTSIG if this is not desired.

• If the ∗Orientation option is used in conjunction with UMAT, stress and strain components will be in the local system (again, this basis system rotates with the material in finite-strain analysis).

• Tensor state variables must be rotated in the subroutine (use ROTSIG).

Usage Hints

• At the start of a new increment, the strain increment is extrapolated from the previous increment.

  • This extrapolation, which may sometimes cause trouble, is turned off with ∗STEP, EXTRAPOLATION=NO

• If the strain increment is too large, the variable PNEWDT can be used to suggest a reduced time increment.

  • The code will abandon the current time increment in favor of a time increment given by PNEWDT*DTIME.

• The characteristic element length can be used to define softening behavior based on fracture energy concepts.

Matmodlab UMAT Support

Matmodlab has built in support for compiling and running Abaqus UMATs

• Coding standards follow that of Abaqus
• All Abaqus utility procedures are implemented
• Develop/Test/Fit the material in Matmodlab before transitioning to the full finite element code

Example: Isotropic Isothermal Linear Elasticity

• Constitutive law:

  $$\pmb{\sigma} = 3K{\rm iso}\pmb{\epsilon} + 2G{\rm dev}\pmb{\epsilon}$$

• Jaumann (corotational) rate form:

  $$\dot{\pmb{\sigma}} = 3K{\rm iso}\dot{\pmb{\epsilon}} + 2G{\rm dev}\dot{\pmb{\epsilon}}$$

• The Jaumann rate equation is integrated in a corotational framework:

The appropriate coding is shown in the following cells.

UMAT Procedure Declaration

subroutine umat(stress, statev, ddsdde, sse, spd, scd, rpl, &
     ddsddt, drplde, drpldt, stran, dstran, time, dtime, temp, dtemp, &
     predef, dpred, cmname, ndi, nshr, ntens, nstatv, props, nprops, &
     coords, drot, pnewdt, celent, dfgrd0, dfgrd1, noel, npt, layer, &
     kspt, kstep, kinc)

  implicit none
  character*8, intent(in) :: cmname
  integer, intent(in) :: ndi, nshr, ntens, nstatv, nprops
  integer, intent(in) :: noel, npt, layer, kspt, kstep, kinc
  real(8), intent(in) :: sse, spd, scd, rpl, drpldt, time, dtime, temp, dtemp
  real(8), intent(in) :: pnewdt, celent
  real(8), intent(inout) :: stress(ntens), statev(nstatv), ddsdde(ntens, ntens)
  real(8), intent(inout) :: ddsddt(ntens), drplde(ntens)
  real(8), intent(in) :: stran(ntens), dstran(ntens)
  real(8), intent(in) :: predef(1), dpred(1), props(nprops), coords(3)
  real(8), intent(in) :: drot(3, 3), dfgrd0(3, 3), dfgrd1(3, 3)

  integer :: i, j
  real(8) :: K, K3, G, G2, Lam
  character*120 :: msg
  character*8 :: charv(1)
  integer :: intv(1)
  real(8) :: realv(1)
  ! ------------------------------------------------------------------------- !

• Free format Fortran and implicit none used
• In Fortran, all variable declarations must appear before any "executabe statements"

Error Checking

if (ndi /= 3) then
     msg = 'this umat may only be used for elements &
          &with three direct stress components'
     call stdb_abqerr(-3, msg, intv, realv, charv)
  end if

• This UMAT is only formulated for plane strain and 3D elements
• Generally speaking, plane stress elements require a different formulation

Stiffness and Stress Update

! elastic properties
  K = props(1)
  K3 = 3. * K
  G = props(2)
  G2 = 2. * G
  Lam = (K3 - G2) / 3.

  ! elastic stiffness
  ddsdde = 0.
  do i=1,ndi
     do j = 1,ndi
        ddsdde(j,i) = Lam
     end do
     ddsdde(i,i) = G2 + Lam
  end do
  do i=ndi+1,ntens
     ddsdde(i,i) = G
  end do

  ! stress update
  stress = stress + matmul(ddsdde, dstran)

  return
end subroutine umat

• Stiffness corresponds to storage of engineering strain components
• Stress is updated in an incremental formulation

The Complete Subroutine

The complete subroutine is written to user_elastic.f90 below

Compile and Run the Model in Matmodlab

• UMATs can be run directly in Matmodlab
• See UserMaterials.ipynb for options for running user materials in Matmodlab

Note The elastic material model is implemented in free-form Fortran. Matmodlab handles free-form Fortran seemlessly. In Abaqus, the user must inform the Abaqus build system that the material model is implemented in free-form Fortran by adding the -free (/free on Windows) flag to the compile_fortran command in their abaqus_v6.env environment file.

• Use the build-fext executable to build the material model. The arguments for build-fext are:

  build-fext <name> <source files>
  
  

  if name is umat, build-fext will compile the subroutine against a library providing several Abaqus subroutines (sprind, std_abqerr, etc).

Model Verification: Uniaxial Stress Step

• Apply axial strain while holding lateral stresses constant
• For a step of uniaxial stress, the slope of stress vs. strain should be the Young's modulus $E$, given by $E=\frac{9KG}{3K+G}$

Model Verification: Uniaxial Strain Step

• Apply axial strain while holding lateral strains constant
• For a step of uniaxial stress, the slope of stress vs. strain should be the constrained modulus $H$, given by $H=K+\frac{4}{3}G$

Conclusion

• Commercial finite element codes' user material interfaces allow custom constitutive responses in simulations
• Most user material interfaces are (still) implemented in Fortran
• Matmodlab natively supports the Abaqus user material interface